Global structure of Bose-Einstein condensates at high rotation: 
Beyond the lowest Landau level description 
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We study the global density profile of a rapidly rotating Bose-Einstein condensate in a harmonic 
trap with transverse frequency ui±. By introducing an additional variational degree of freedom 
to the lowest Landau level wave function, we demonstrate that with increasing strength of the 
interparticle interaction, the global density profile changes from a Gaussian to the inverted parabolic 
one characteristic of Thomas-Fermi theory. The criterion for the lowest Landau level wave function 
to be a good approximation for the global structure is that the mean field energy be small compared 
with hcj±/N v , where N v is the number of vortices in the cloud. This condition is more stringent 
than the requirement that the mean field energy be small compared with hui± which is necessary 
for the lowest Landau level wave function to be a good approximation to the local structure. Our 
results show that the lowest Landau level wave function is inappropriate for the global structure 
of the system realized in recent experiments even though this wave function can describe the local 
structure well. 

PACS numbers: 03.75.Hh, 05.30.Jp, 67.40.Vs, 67.40.Db 
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Quantized vorticity and vortex lines are characteristic 
features of superfluids Q, @ and since the first observa- 
tion of a vortex in a Bose-Einstein condensed atomic gas 
many beautiful experiments have been performed to 
observe vortex lines and vortex lattices in these systems 
000 One fundamental question is how condensates 
behave when the vortex core size become comparable to 
the spacing between vortices. Such conditions have not 
been achieved for superfluid liquid helium 4, but they 
have been for atomic gases in harmonic traps. Ho noted 
that the Hamiltonian for a rotating gas in a harmonic 
trap is similar to that for charged particles in a magnetic 
field, and he argued that for rotational angular veloci- 
ties just below the radial trap frequency uj±_ all particles 
would condense into the lowest Landau level (LLL) of 
the rotational motion 0. Motivated by this suggestion, 
Schweikhard et at have recently achieved rotational an- 
gular velocities f2 in excess of 0.99wj_, at which the cloud 
contains a number of vortices of order one hundred |lfj| . 
The frontiers of experiments have reached the so-called 
"mean-field quantum Hall" regime where the Ml is large 
compared with the interaction energy, that is Ml ng, 
where n is the particle density and g = AnTi 2 a/m is the ef- 
fective two-body interaction, m being the particles mass 
and a the scattering length. 

Employing a quantum-Hall like wave function, which 
contains only components in the lowest Landau levels, Ho 
predicted that the smoothed density profile of a trapped 
cloud would be Gaussian. To understand the difference 
between the conventional approach to slowly-rotating 
condensates and rapidly-rotating ones, the authors of 
Ref. (see also Ref. 01) adopted a more general wave 
function consisting of the product of a slowly varying en- 
velope function, which describes the global structure of 
the cloud and a rapidly varying one which describes the 



properties of individual vortices. They concluded that 
the density profile would have the inverted parabolic form 
characteristic of Thomas-Fermi theory, rather than the 
Gaussian one. 

In the present paper, we consider the properties of a 
rotating cloud in a potential V(r) = muij_r 2 /2, using a 
trial wave function which is a generalization of the lowest 
Landau level one. We find that the global structure of 
the cloud is of the Thomas-Fermi form provided Ml 3> 
ng/N v , where N v is the number of vortices within the 
cloud. For simplicity, we confine ourselves mainly to the 
two-dimensional problem in this paper. 

We shall use for the wave function 'F of the condensed 
state the expression 



<F(r) = JV^tO&xl = iV 1/2 V(r) , 



(1) 



where N is the number of particles in the condensate 
and the lowest Landau level wave function </>lll has the 
general form (i.e., quantum-Hall like wave function) as 
in Ref. 
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where £ = x + iy, d are the vortex positions, a± = 
(H/mLdjJ 1 / 2 and is the normalization constant. The 
function h which changes the density profile from the 
form predicted by the lowest Landau level calculation, 
is real, and varies slowly on the scale of the intervortex 
separation. In the lowest Landau level wave function, 
the positions of vortices determine the density distribu- 
tion, apart from an overall multiplicative constant, and 
the introduction of the modulating function allows us to 
break this requirement. The normalization condition for 



2 



the wave function is J <i 2 r/i 2 |</>LLL| 2 = 1- The normaliza- 
tion condition for <^>lll is arbitrary, but for dcfinitcncss 
we shall make the choice J d 2 r|</>LLL| 2 = 1. 

The energy per particle of the condensate in the ro- 
tating frame is given by E' = E — CIL Z , where E is the 
energy in the non-rotating frame and L z is the expecta- 
tion value of the angular momentum per particle about 
the axis of rotation. Following Ho |9(, we write this in 
the form 



E' = (ljx - tyL z 



+ / d z r ib* 



a;j_x r 
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where wi = wj_z. 

An important observation is that even if the spatial 
variations of h are of order unity over the cloud, the ad- 
mixture of excited Landau levels in the wave function 
JTJ is of order a±dh/dr ~ a±/R ~ relative to 

the LLL contribution [Tjj. Here R is the radial extent 
of the cloud and N v ~ R 2 ja\ is the number of vortices 
in the cloud. From this one can show that, apart from 
corrections of order h(tu± — Q) /N v , for the wave function 
(JTJ, one finds 



E' = hfl 



d 2 r 



LLL 



+ ^±-tt)^ 2 + ^(|0LLL|> 4 
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Here <?2D is the effective coupling parameter in two di- 
mensions and b = ( | </>lll | 4 ) / ( | </>lll | 2 ) 2 is a factor of order 
unity describing the renormalization of the effective in- 
teraction due to the rapid density variations on the scale 
of the vortex separation 

mill 

If the wave function is uniform in the direction of the 
axis of rotation, <?2D = Ng/Z, where Z is the axial ex- 
tent of the cloud, while if the wave function for motion 
in the direction of the rotation axis corresponds to the 
ground state of a particle in a harmonic oscillator po- 
tential with frequency uj z , g 2 D — Ng/(\f2~jra z ), where 
fl 2 = (h/mwz) 1 / 2 Since the terms in braces in Eq. 

flUl vary slowly in space, we may obtain a good approx- 
imation to the energy by making the averaged vortex 
approximation as in Ref. 0. Like Ho, we shall in this 
paper assume that the number of vortices per unit area 
is uniform. If the vortex density is nonuniform, this will 
increase the interaction energy, which for a uniform sys- 
tem is a minimum for a regular triangular lattice. For 
such a uniform array of vortices, one finds 
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Here (. . .) denotes an average over an area of linear size 
large compared with the vortex separation but small 
compared with a. The width a is given by 1/er 2 = 



l/<2j_ — nn v , where n v is the average vortex density in 
the plane perpendicular to the rotation axis [tj. 

The first term in braces in Eq. represents the ex- 
tra kinetic energy associated with admixture of excited 
Landau levels. It scales with the radius of the system in 
the same way as the interaction energy. If g 2 T> *C h 2 /m, 
the first term suppresses spatial variations of h, and the 
solutions is h « 1, and the wave function reduces to the 
LLL one. The condition for the LLL wave function to be 
a good approximation may also be written as Na/Z <C 1. 
In current experiments Na/Z is of order 10-100, so this 
condition is strongly violated. In the opposite limit, 
52D ^ h 2 /m, the extra kinetic energy is unimportant, 
and the optimal density profile is obtained by minimiz- 
ing the second and third terms in the integrand, which 
results in a Thomas-Fermi density profile. Another way 
to express the criterion for validity of the LLL approxi- 
mation is that the interaction energy per particle, gn, be 
small compared with hiv±/N v , rather than the condition 
Tiu!± which has generally been assumed in earlier work. 

To obtain quantitative results we perform a variational 
calculation with a trial function that interpolates be- 
tween the Gaussian and Thomas-Fermi forms. We ex- 
ploit the fact that lim a ^ oc (l — t/a) a = e~ l as Fetter did 
in calculations for non-rotating clouds 0] . Thus we take 



h(r) =A h [l 



2 /2a 2 



(6) 



for < r < ^J~olL and otherwise h = 0. The normaliza- 
tion factor A h is given by A 2 h = (er 2 /L 2 )(l + 1/a). The 
number of particles per unit area, divided by TV, is given 



by v{r) = h 2 ^: 



lll| 2 ) with the above expression (JSJ for 
h can describe both Gaussian form (a — * 00) and the 
Thomas- Fermi one (a = 1). Hereafter, we shall refer to 
the expression @) for h as the Gaussian - Thomas-Fermi 
(G-TF) form. 

Substituting Eq. © into Eq. , we obtain the energy 
Eq_ tf for the Gaussian-TF modulation 
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Minimizing E' G _ TF with respect to a and L, we obtain 
the following expressions for these two quantities: 
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is a dimensionless parameter determining the strength of 
interparticle interactions. 

Then we minimize E' G _ TF with respect to a. The op- 
timal value of a is determined by the quartic equation 

I 



a 



since the three other solutions give negative or complex 
values of a. In the limits of weak (k — * 0) and strong 
interaction (k — * oo) the above solution behaves as 

a ~ 8/k (k < 1) , and (12) 

3 1 2 1 

a ~ 1 + («^»1). (13) 

Observe that a does not depend on f2. 

The shape index a given by Eq. (|llfl is plotted in Fig.^ 
and it decreases from infinity in the absence of interaction 
to unity as the strength of the interaction (k) increases 
from zero to infinity. This means that the density pro- 
file of the cloud changes from a Gaussian to an inverted 
parabola due to the interaction. Without the minimiza- 
tion with respect to a, we can describe the cloud with 
an assumption on the shape of the density profile corre- 
sponding to the given value of the shape index. For the 
parameters L and <r, we can use the expressions of Eqs. 
© and 0- If we take a — ► oo, we reproduce a Gaus- 
sian density profile as in the case without the modulating 
function, but its width is optimized. 

Now we compare the results of the following two cases: 
the Gaussian-TF profile with an optimized value of a 
and the Gaussian profile with a — > oo, and discuss the 
situation relevant to the experiments. For definiteness, 
we write the energy for the Gaussian profile as E' G , which 
can be obtained by taking a — > oo in Eq. Q. Relative 
energy differences (E' G _ TF — E' G )/(E' G — hfl) between the 
Gaussian-TF and Gaussian cases are shown in Fig. [3 In 
the limit k — > oo, the leading terms are 

i 

e' g ~ nn + Vz huj^ (i - 2 kS , (14) 

4 / , 

e' g TF ~ nn + -huj^ i k3 . (15) 

3 \ wi/ 

Thus, in this limit, the relative energy difference 
(E' G _ TF -E' G )/(E' G - nSl) converges to (2 3 / 2 - 3)/3 ~ 
—0.0572 irrespective of the value of f2 (see Fig. |2J. 
This result shows us that the density profile can change 



a 4 - 2(1 + 4A)a 3 - 12Aa 2 + 2(1 - 3A)a - 1 - A = 0, where 
A = K . Solving this equation for a, we find that the 
only physical solution is 



(11) 
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FIG. 1: Shape index a as a function of kT 1 covering a range 
of interaction strengths typical of current experiments (k ~ 
100). In the inset, a is plotted over a wider range of to 
show its asymptotic behavior 1121 . 



markedly even though the energy reduction is less than 
6%. 

We now consider the system realized in Ref. 10] rotat- 
ing at fl ~ 0.99w^, which is typical of the highest angular 
velocities achieved so far. If we assume that the density is 
uniform in the axial direction, we get K = (2Na)/Z > 52 
using the relation <?2D = gN/Z and the experimental data 
Na/Z > 26 given in Ref. [T0| - If we assume a Gaussian 
density profile in this direction instead of the uniform 
one, we obtain k — 2Na/ (\/2Tra z ) ~ 1.4 x 10 2 using 
92B = gN/(V2~Tra z ), oj z = 2ir x 5.3 Hz, N = 1.5 X 10 5 
at — 0.989o;^, and a — 5.6 nm for the triplet state 
of 87 Rb. Thus we conclude that the typical value of the 
parameter k relevant to the experiments of Ref. |Toj is of 
order 10 2 . 

In Fig. |3J we plot the global profile for the number of 
particles per unit area v(r) — /i 2 (|^lll| 2 ) at n = 0.99wj_ 
and k — 100 for the Gaussian-TF case in addition to 
those for the two extreme cases of the Gaussian form with 



\{l + 4A) + 1^/(1 + 4A) 2 + 8A - 3 ■ 2§ [A(l + A)]* 
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FIG. 2: Difference between the energies E G _ TF calculated in 
the Gaussian-TF case with the optimized a given by Eq. If ft 
and E' G in the Gaussian case with a — > oo, relative to E' G — hQ. 
For any value of f2, the curve converges to (2 3 ^ 2 — 3)/3 ~ 
—0.0572 in the limit of k — > oo. 
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FIG. 3: Global density profile v(r) = h 2 (\ <^lll| ) of the cloud 
for £1 — 0.99u± and k = 100. The solid line is obtained for the 
optimized a given by Eq. If ft . The dotted line and dashed 
one are obtained for a — > oo and a = I , respectively. 



a — > oo and the inverted parabolic one labelled as "TF" . 
The latter is obtained by the procedure corresponding 
to the Thomas-Fermi approximation in the non-rotating 
case in which we neglect the contribution from the first 
term in the braces in Eq. [i.e., the second term in Eq. 
(JZJ)]. As can be seen from the solid line of this figure, the 
density profile for the Gaussian-TF modulation is closer 
to the TF case than to the Gaussian one and it can be 
fitted rather well b y a n inverted parabola except at the 
edge of the cloud 17]. From Fig. one can see that 
the relative energy reduction compared with that for the 
Gaussian approximation is only ~ 2.9%. 

In this paper, we have demonstrated that the global 
density profile of a rapidly rotating Bose-Einstein con- 
densate is not well described by the lowest Landau level 
wave function, even though the mean field energy is small 
compared with huj±, a condition satisfied in recent ex- 
periments. By performing a variational calculation with 
a trial wave function which includes components in ex- 
cited Landau levels, we have shown how the density pro- 
file changes from a Gaussian to the Thomas-Fermi form 
as the strength of the interparticle interaction increases. 
The reduction in the energy resulting from the extra de- 
gree of freedom in the wave function is of the order of 
one per cent. The approach described in this paper may 
be generalized to time-dependent problems by deriving 
equations for the evolution of the variational parameters 
from the condition that the action be stationary. 

This work arose as a result of discussions and corre- 
spondence with Gordon Baym, to whom we are very 
grateful. In addition, we thank V. Cheianov for help- 
ful comments. G. W. thanks K. Sato, K. Yasuoka, and 
T. Ebisuzaki for continuous encouragement and L. M. 
Jensen, and P. Urkedal for arranging the computer en- 
vironment. This work was supported by Grants- in- Aid 
for Scientific Research provided by the Ministry of Edu- 
cation, Culture, Sports, Science and Technology through 
Research Grant No. 14-7939. 



[1] R. J. Donnelly, Quantized Vortices in Helium II (Cam- 
bridge University Press, Cambridge, 1991). 

[2] E. B. Sonin, Rev. Mod. Phys. 59, 87 (1987). 

[3] M. R. Matthews, B. P. Anderson, P. C. Haljan, D. S. 
Hall, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. 
83, 2498 (1999). 

[4] A. L. Fetter and A. A. Svidzinsky, J. Phys. : Condens. 
Matter 13, R135 (2001). 

[5] K. W. Madison, F. Chevy, W. Wohlleben, and J. Dal- 
ibard, Phys. Rev. Lett. 84, 806 (2000). 

[6] J. R. Abo-Shaeer, C. Raman, J. M. Vogels, and W. Ket- 
terle, Science 292, 476 (2001). 

[7] P. C. Haljan, I. Coddington, P. Engels, and E. A. Cornell, 
Phys. Rev. Lett. 87, 210403 (2001). 

[8] P. Engels, I. Coddington, P. C. Haljan, and E. A. Cornell, 
Phys. Rev. Lett. 89, 100403 (2002). 

[9] T.-L. Ho, Phys. Rev. Lett. 87, 060403 (2001). 



[10] V. Schweikhard, I. Coddington, P. Engels, V. P. Mogen- 
dorff, and E. A. Cornell, cond-mat/0308582 

[11] G. Baym and C. J. Pethick, Phys. Rev. A (in press), 
(cond-mat/0308325). 

[12] U. R. Fischer and G. Baym, Phys. Rev. Lett. 90, 140402 
(2003). 

[13] This will be demonstrated in detail elsewhere (G. Watan- 

abe and C. J. Pethick, in preparation). 
[14] J. Sinova, C. B. Hanna, and A. H. MacDonald, Phys. 

Rev. Lett. 89, 030403 (2002). 
[15] A. L. Fetter, Phys. Rev. A 68, 063617 (2003). 
[16] A. L. Fetter, J. Low Temp. Phys. 106, 643 (1997). 
[17] In the region r < 10ax , the solid curve in Fig. [3] can 

be fitted quite well by the inverted parabola v{r)a\ — 

0.00343 - 2.1 x W- 5 (r/a 1 _) 2 . 



